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ABSTRACT 


In  satellite  tracking  using  ground-based  radars,  an  estimate  of  the  total  electron  content  (TEC) 
along  the  path  to  the  satellite  is  required  to  measure  accurately  the  range  of  the  satellite.  This 
estimate  is  necessary  because  the  radar  wave  travels  at  a  slower  speed  as  it  propagates  through  the 
ionosphere.  The  range  error  A R  that  is  introduced  is  dependent  on  the  radar  frequency /and  on  the 
TEC  along  the  propagation  path,  and  can  be  expressed  by 
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where  N  is  the  local  electron  density  and  R  is  the  radar  range.  The  TEC  can  vary  significantly  with 
the  time  of  day,  geomagnetic  activity,  and  look  direction.  A  real-time  synoptic  ionospheric 
monitoring  system  has  been  developed  using  data  acquired  from  a  TI4100  Global  Positioning 
System  (GPS)  receiver  for  use  at  the  Millstone  Hill  satellite-tracking  radar.  The  TI4100  GPS 
receiver  can  track  up  to  four  GPS  satellites  at  any  one  time.  Each  GPS  satellite  transmits  signals  at 
two  different  L-band  frequencies:  LI  (1575.42  MHz)  and  L2  (1227.6  MHz).  The  TEC  along  the 
path  to  each  satellite  can  be  determined  by  combining  both  frequencies  using  the  pseudorange  and 
the  integrated  phase  data.  At  Millstone,  the  TEC  is  measured  every  3  s  for  each  GPS  satellite  in  view. 
The  data  are  input  into  a  Kalman  filter  that  is  used  to  predict  the  coefficients  of  a  simple  TEC  model 
with  azimuth  and  elevation  dependence.  This  model  takes  advantage  of  the  real-time  knowledge 
provided  by  the  GPS  data  of  the  variations  in  TEC  around  the  Millstone  location.  The  coefficients 
for  this  model  are  then  sent  to  the  satellite-tracking  computer,  and  the  model  is  applied  in  real  time 
to  account  for  the  ionospheric  path  delay  to  whatever  satellite  is  currently  in  track.  The  preliminary 
results  of  using  this  ionospheric  monitoring  system  at  Millstone  will  be  discussed.  The  zenith  value 
of  the  TEC  predicted  by  our  GPS  model  will  be  compared  with  another  ionospheric  measurement, 
the  foF2  obtained  from  the  collocated  University  of  Lowell  Digisonde.  The  TEC  values  predicted 
by  our  GPS  model  during  both  geomagnetically  quiet  and  disturbed  time  periods  will  be  discussed, 
as  well  as  those  associated  with  high  solar  flux  time  periods.  Finally,  the  improvement  in  our  radar 
system  calibration  due  to  the  use  of  the  new  GPS  model  will  be  demonstrated.  This  improvement 
is  evident  in  observations  of  the  average  range  residual  over  a  pass  of  the  Lageos  satellite.  The 
standard  deviation  of  these  residuals  drops  from  approximately  20  TEC  units  to  about  5  TEC  units 
when  the  new  GPS  model  is  used  to  estimate  the  ionospheric  refraction  correction.  A  TEC  unit  is 
1016  electrons/m2.  At  the  Millstone  radar  frequency  of  1295  MHz,  this  correction  corresponds  to 
an  improvement  in  the  standard  deviation  from  5.3  to  1.4  m.  It  is  our  opinion  that,  other  than 
upgrading  all  satellite-tracking  radars  to  dual-frequency  capability,  the  GPS  is  currently  the  best 

monitoring  system  of  the  TEC  available.  _ _ 
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1.  INTRODUCTION 


A  transportable  device  has  been  developed  at  Lincoln  Laboratory  that  uses  the  Global 
Positioning  System  (GPS)  to  determine  real-time  ionospheric  path  delays  for  use  in  satellite 
tracking.  This  system  was  integrated  into  the  Millstone  Hill  real-time  satellite-tracking  software  in 
March  1991.  This  report  will  present  the  initial  satellite-tracking  results  using  this  system  and  will 
describe  the  Kalman  filter  used  in  obtaining  the  model  coefficients.  Before  discussing  the  real-time 
data  processing  involved,  some  background  information  will  be  presented  on  both  the  ionospheric 
path  delay  and  on  the  design  and  use  of  the  GPS.  This  material  will  illustrate  the  importance  of 
correctly  estimating  the  ionospheric  path  delay  and  will  explain  how  the  GPS  data  can  be  used  to 
determine  it. 


The  ionospheric  path  delay  is  a  consequence  of  the  slower  group  velocity  v^r  of  the  radar  wave 
as  it  travels  through  the  ionosphere.  Therefore,  the  true  range  to  the  satellite  Rr  is  somewhat  less 
than  the  apparent  range  Ru,  which  assumed  that  the  radar  wave  traveled  at  the  speed  of  light  c  in  a 
vacuum. 

A R  =  R-R=ct-vt.  (2) 

a  t  gr 


The  error  in  the  measurement  of  radar  range  A R  is  dependent  on  the  radar  frequency/and  the 
total  electron  content  (TEC)  along  the  ray  path  from  the  radar  to  the  satellite.  It  can  be  expressed 
by  the  following  equation: 


A R  (meters)  = 


40.3 

f2 


(3) 


where  Ne  is  the  local  electron  density,  and  the  integrated  value  of  it  is  the  TEC.  The  TEC  can  vary 
with  look  angle,  time  of  day,  season,  geomagnetic  activity,  and  solar  flux  value.  Typical  errors  at 
L-band  frequencies  in  radar  range  due  to  the  ionosphere  are  about  32  m  at  20  deg  of  elevation.  At 
UHF  frequencies,  errors  are  estimated  to  be  280  m  at  20  deg  of  elevation,  nearly  1 0  times  larger  than 
at  L-band. 


The  GPS  design  included  the  broadcast  of  two  different  frequencies  so  that  the  ionospheric 
delay  term  could  be  eliminated.  Each  GPS  satellite  is  broadcasting  a  signal  at  LI  (1575.42  MHz) 
and  at  L2  (1227.6  MHz).  By  computing  the  difference  in  delay  of  the  LI  and  L2  signals,  the 
ionospheric  term  can  be  computed  directly.  In  practice,  there  is  also  the  additional  problem  of 
multipath. 

Unlike  most  GPS  users,  we  do  not  consider  the  ionospheric  delay  term  to  be  noise.  Rather, 
we  are  using  the  GPS  to  measure  the  ionospheric  TEC  along  the  path  to  each  GPS  satellite  in  view. 
The  Millstone  GPS  Real-Time  Ionospheric  Mapping  System  (GRIMS)  [1]  monitors  the  local 
ionosphere  by  continuously  sampling  both  L-band  frequencies  for  each  satellite  in  view.  The  group 
delay  and  reconstructed  carrier  waves  are  then  used  to  construct  the  TEC.  By  combining  the  TEC 
measurements  from  all  satellites,  a  four-dimensional  model  (space  and  time)  of  the  ionosphere  is 
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generated.  The  coefficients  of  this  model  are  then  sent  to  our  satellite-tracking  computer  and  are 
used  in  real  time  to  compute  the  ionospheric  correction  in  the  direction  that  the  Millstone  radar  is 
pointing. 

All  results  presented  here  are  given  in  TEC  units,  the  number  of  electrons  x  1 0,6/m2.  For  easy 
reference,  1  m  of  range  correction  at  LI  =  6.15929  TEC  units.  1  TEC  unit  =  2.85  x  TDD,  where 
TDD  is  equal  to  TTL2  -TTL 1  (transit  time  of  L2  -  transit  time  of  L 1  in  ns),  and  stands  for  time  of 
differential  delay. 
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2.  GPS  TEC  MEASUREMENTS 


In  principle,  the  ionospheric  path  delay  can  be  obtained  directly  from  the  difference  in  the  GPS 
pseudorange  measurements.  However,  although  the  GPS  group  delay  measurements  are  accurate 
enough  for  the  majority  of  user  needs,  the  ionospheric  correction  cannot  be  determined  directly  from 
this  data.  This  effect  is  due  primarily  to  multipath  and  the  L1/L2  biases. 

Multipath  effects  in  the  group  delay  data  account  for  approximately  a  6-ns  increase  in  the  noise 
level.  For  comparison,  a  2-ns  increase  in  group  delay  amounts  to  a  1-m  error  in  the  path  delay  at 
L-band.  At  Millstone,  the  multipath  effects  have  been  reduced  by  the  use  of  absorbing  material 
under  the  GPS  antenna  and  can  be  further  reduced  by  data  averaging  [2].  The  multipath  problem 
is  not  as  severe  in  the  GPS  integrated  phase  measurements.  The  GPS  phase  data  show  an  increase 
in  the  noise  level  on  the  order  of  only  50  ps  [3].  The  trouble  with  the  phase  data  is  that  it  contains 
an  unknown  bias  due  to  the  integer  cycle  ambiguity  of  phase  measurements  (i.e.,  the  number  of  full 
phase  cycles  between  the  satellite  and  receiver  is  not  known).  By  combining  the  GPS  integrated 
phase  data  with  the  group  delay  data,  the  difficulties  inherent  with  both  data  types  may  be  overcome. 

Figure  1  illustrates  both  the  GPS  integrated  phase  and  GPS  group  delay  data  plotted  as  a 
function  of  time  for  one  satellite  (SSC  number  1 1054,  PRN  number  6)  that  was  in  view  from 
Millstone  during  a  3-h  time  span  in  March  1989.  The  increased  noise  level  due  to  multipath  is 
apparent  in  the  group  delay  data.  The  smooth  appearance  of  the  integrated  phase  data  shows  that 
it  is  less  susceptible  to  multipath.  To  account  for  the  unknown  cycle  ambiguity  and  bias,  the  mean 
value  of  the  phase  data  is  set  equal  to  the  mean  value  of  the  group  delay  data. 


Figure  I.  Illustration  of  GPS  phase  and  group  delay  data.  Satellite  SSC  no.  11054,  PRN  no.  6.  /  March  1989. 
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Another  problem  with  making  TEC  measurements  using  GPS  data  is  that  both  the  satellite  and 
the  receiver  hardware  introduce  an  unknown  bias — or  additional  delay — between  the  LI  and  L2 
signals.  These  biases  must  be  measured  or  estimated  and  removed  from  the  GPS  group  delay  data 
in  order  to  determine  accurately  the  ionospheric  path  delay.  These  biases  have  been  studied 
extensively  [3-6]. 

Coco  et  al.  [5]  have  shown  that  although  the  TI4100  receiver  (the  GPS  receiver  used  at 
Millstone)  bias  varies  from  unit  to  unit,  the  change  over  time  in  the  bias  is  quite  small.  They 
estimated  the  day-to-day  variation  of  the  TI4 1 00  bias  over  a  five-week  period  to  be  less  than  0.5  ns 
(corresponding  to  a  23-cm  effect  at  LI  or  1.4  TEC  units).  The  Millstone  TI4100  GPS  receiver  is 
assumed  to  have  a  negligible  bias  because,  at  times,  the  difference  between  the  Millstone  GPS  TEC 
measurements  and  those  of  the  collocated  incoherent  scatter  (IS)  radar  has  been  measured  to  be  less 
than  1  TEC  unit  (0.35  ns  of  differential  delay). 

Satellite  biases  are  more  difficult  to  determine  and  pose  a  serious  problem  to  ionospheric 
estimation  when  using  the  GPS.  Studies  have  developed  techniques  for  determining  the  space 
vehicle  (SV)  biases  using  a  multistation  system  of  receivers  [3,4],  The  JPL-determined  biases  [7] 
have  been  adopted  for  the  work  reported  here,  as  this  is  the  standard  convention  in  the  ionospheric 
community. 

In  Figure  2,  a  comparison  of  GPS,  incoherent  scatter,  and  Hamilton  TEC  measurements  is 
shown.  The  measurements  were  taken  on  27  April  1990.  The  GPS  receiver  at  Millstone  was 
tracking  the  GPS  satellite  PRN  number  6  while  it  was  in  view.  At  the  same  time,  the  Millstone  UHF 
IS  radar  was  pointed  so  that  it  was  taking  IS  measurements  of  the  ionosphere  along  the  same  line 
of  sight  to  this  GPS  satellite.  The  IS  data  were  then  processed  to  give  an  alternate  determination  of 
the  TEC  [8].  The  IS  data  were  carefully  calibrated  at  the  beginning  and  end  of  each  run  using  foF2 
data  from  a  collocated  Digisonde  that  was  operated  by  the  University  of  Lowell.  The  Millstone  UHF 
radar  can  provide  electron  density  profiles  from  approximately  100  to  1000  km  or  higher.  Due  to 
the  decreasing  signal-to-noise  ratio,  however,  data  at  altitudes  above  800  km  have  been  eliminated 
here.  The  IS  TEC  measurements  presented  in  this  plot  have  been  derived  by  integrating  the  electron 
density  profiles  from  100  to  800  km  in  height. 

Figure  2  aLo  shows  a  comparison  of  our  GPS  data  with  the  TEC  data  taken  at  Hamilton, 
Massachusetts.  Using  the  technique  of  Faraday  rotation,  the  polarimeter  at  Hamilton  measures  the 
change  in  phase  of  a  signal  transmitted  by  the  GOES-B  satellite.  The  Faraday  rotation  technique 
measures  the  TEC  up  to  approximately  2500  km  [7],  due  to  the  smaller  magnetic  field  at  higher 
altitudes.  GOES-B  is  at  a  range  of  approximately  39,000  km.  Because  Hamilton  is  at  nearly  the  same 
latitude  and  at  a  comparable  longitude  to  Millstone,  our  GPS  measurements  of  TEC  can  be  compared 
to  Hamilton’s  measurements  as  long  as  our  measurements  are  made  looking  along  a  similar  line  of 
sight  (LOS)  to  that  of  Hamilton's  LOS  to  GOES-B.  The  comparison  shown  in  Figure  2  has  made 
some  adjustment  for  the  different  elevations  of  the  GPS  and  Hamilton  measurements.  The 
difference  between  the  Hamilton  and  GPS  TEC  data  is  less  than  2  TEC  units.  Basically,  Figure  2 
compares  the  GPS  measurement  of  TEC  out  to  19,000  km,  the  Faraday  rotation  measurement  out 
to  2500  km,  and  the  IS  measurement  of  TEC  out  to  800  km.  The  large  difference  between  the  IS 
measurements  of  TEC  and  the  Hamilton  and  the  GPS  measurements  of  TEC  can  be  attributed  to  me 
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fact  that  the  ionosphere  extends  significantly  above  800  km  during  time  periods  of  high  solar  flux. 
In  April  1990,  the  FI 0.7  cm  flux  was  more  than  200.  Further  information  on  TEC  observations 
detected  above  800  km  can  be  found  in  associated  literature  [9]. 


201711-2 


Figure  2.  Comparison  of  GPS,  incoherent  scatter,  and  Hamilton  TEC  measurements.  GPS  Satellite  SSC  no.  I  I054.PRN 
no.  6,  27  April  1990. 
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3.  GPS  PREPROCESSING 


The  GPS  data  is  preprocessed  before  it  is  sent  to  the  Kalman  filter.  This  preprocessing  includes 
reconstructing  the  LI  and  L2  SV  times,  applying  the  SV  cioc  k  corrections,  and  fixing  detected  cycle 
slips  in  the  integrated  phase  data.  The  algorithms  used  in  the  preprocessing  are  described  here. 

The  TI4100GPS  receiver  acquires  the  GPS  data.  Then,  in  the  preprocessing,  the  raw  LI  and 
L2  SV  times  are  reconstructed,  and  the  space  vehicle  (SV)  clock  corrections  are  applied.  The  SV 
clock  correction  is  given  by 

2 

At  =  afr.+af,(t-t  )  +  afJt-t  )  +A/  ,  (4) 

sv  /  0  /1\  oc)  f2\  oc)  r 

where 


A /  =  the  satellite  (SV)  clock  correction  (s), 

ay0  =  polynomial  coefficient  (s), 
djy  =  polynomial  coefficient  (s/s), 

=  polynomial  coefficient  (s/s2), 
t  =  GPS  system  time  (s), 
toc  =  clock  data  reference  time  (s), 

A tr  -  relativistic  correction  term  (s). 

In  addition. 


where 


(5) 


R  =  instantaneous  position  vector  of  the  SV  (m), 

V  =  instantaneous  velocity  vector  of  the  SV  (m/s), 
c  =  speed  of  light  (m/s). 

The  navigation  (NAV)  data  are  used  to  calculate  the  azimuth  and  elevation  for  each  data  point. 

Sometime  '  the  GPS  receiver  loses  lock  on  a  satellite.  When  this  happens,  the  receiver  loses 
track  of  the  phase  front  cycle  count  for  the  satellite  in  question.  This  condition  is  called  a  cycle  slip. 
To  correct  for  cycle  slips,  the  phase  range,  pseudorange,  and  SV  time  data  are  entered  into  a  routine 
that  tests  for  (and,  if  necessary,  fixes)  cycle  slips. 

The  algorithm  used  for  fixing  cycle  slips  is  based  on  the  JPL  GPS  data-editing  algorithm 
described  by  Blewitt  [10].  This  method  has  the  important  advantage  that  it  can  be  implemented  in 
real  time.  The  technique  uses  linear  combinations  of  the  recorded  carrier  phases  and  the  P  code 
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pseudoranges  to  compute  a  wide-lane  bias.  Jumps  in  the  wide-lane  bias  indicate  data  outliers  and 
cycle  slips.  Outliers  are  thrown  out,  and  cycle  slips  are  fixed  by  determining  the  integer  offset 
between  the  data  arcs  before  and  after  the  cycle  slip. 

The  linear  combinations  of  the  phase  and  pseudorange  residuals  were  found  to  be  the  most 
useful  in  our  system.  The  routine  implemented  by  the  T14 1 00  process  uses  these  data  to  determine 
whether  a  cycle  slip  has  occurred.  When  two  consecutive  measurements  for  an  SV  indicate  the 
presence  of  a  cycle  slip,  then  the  algorithm  attempts  to  correct  the  cycle  slip. 

If  the  measurements  are  accepted  by  the  cycle  slip  algorithm,  then  the  azimuth,  elevation,  and 
range  are  computed  from  the  ephemeris.  Range  residuals  are  produced  for  both  the  phase  and 
pseudorange  data. 

The  measurements  that  pass  the  cycle  slip  test  are  then  passed  on  to  the  Kalman  filter  sequential 
estimator  routines.  The  filter  uses  as  input  the  fractional  part  LI  and  L2  SV  times  'in  seconds),  the 
x-y-z  position,  and  the  LI  and  L2  phase  ranges  for  a  single  satellite.  The  Kalman  filter  estimates 
the  current  filter  state  of  the  six  coefficients  that  represent  the  local  ionosphere.  The  model  is 
described  in  the  next  section.  The  Kalman  filter  employs  a  fading  memory  mechanism  that  ensures 
that  the  most  recent  data  are  weighted  more  heavily  than  older  data.  In  addition,  a  second  cycle-slip 
test  is  applied  by  the  Kalman  filter. 
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4.  GPS  REAL-TIME  MODEL 


The  model  of  the  ionosphere  used  in  our  real-time  GPS  system  has  two  components:  a  vertical 
TEC,  which  depends  on  both  azimuth  and  elevation  from  the  site,  and  a  mapping  function,  which 
depends  on  elevation  and  range.  The  mapping  function  translates  the  vertical  TEC  value  into  a  TEC 
value  that  is  valid  along  the  line  of  sight.  This  mapping  function  is  described  in  Appendix  A.  In 
solving  for  the  vertical  TEC  term,  six  coefficients  are  estimated,  a0  through  ay  These  six  coeffi¬ 
cients  are  the  parameters  that  define  the  ionosphere  around  the  Millstone  location.  The  model  for 
the  zenith  correction  has  a  center  point,  or  pole,  defined  at  a  90°  elevation.  The  pole  corresponds 
to  a  zenith  correction  term  aQ  that  is  the  estimate  for  the  ionospheric  TEC  directly  above  Millstone. 
All  of  the  TEC  measurements  are  used  in  calculating  the  zenith  TEC  v  due  over  Millstone.  The  area 
around  the  pole  is  then  divided  in  azimuth  into  quintets.  Separating  each  quintet  is  a  line  that 
represents  the  linear  interpolation  of  the  vertical  TEC  from  an  elevation  equal  to  90°  down  to  0°. 
These  lines  are  defined  by  the  following  equation  that  depends  on  both  the  elevation  and  on  the 
coefficients,  ajt  i  =  0 . 5. 

TEC( vertical )  =  aQ  +  a.  .'  =  1,5  .  (6) 

These  quintets  represent  five  different  azimuth  regions  of  the  ionosphere  around  Millstone,  as 
seen  in  Figure  3. 


2017.1-3 

8-1  =  0° 


Figure  3.  Circus  Tent  model  of  ionosphere  around  Millstone  location. 
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The  center  of  Figure  3  corresponds  to  a  90°  elevation  above  Millstone;  however,  the  rim  of 
the  circle  is  at  a  0°  elevation.  The  pole  aQ  that  is  at  a  90°  elevation,  represents  our  estimate  of  the 
TEC  directly  above  Millstone.  The  zero  azimuth  (a,)  relative  to  this  model  is  offset  from  the 
geographical  0°  azimuth  by  346°  to  align  the  relative  coordinate  system  with  the  geomagnetic  north 
pole.  This  alignment  was  done  because  geomagnetic  coordinates  provide  a  more  natural  system  of 
describing  the  ionosphere.  For  example,  an  electron  density  trough  forms  frequently  at  night  toward 
the  north  from  Millstone,  and  this  trough  tends  to  align  itself  around  the  geomagnetic  north  pole. 


The  advantage  of  this  parameterization  of  the  ionosphere  is  that  it  preserves  the  azimuth  and 
elevation  information  that  is  available  from  the  different  GPS  satellites  in  view.  If  no  GPS  satellite 
is  in  view  in  either  of  the  two  neighboring  quintets,  then  the  aQ  value  is  used  to  estimate  the  TEC 
in  that  direction.  In  the  real-time  satellite-tracking  program,  a  linear  interpolation  in  azimuth 
between  'wo  neighboring  lines  is  used  to  compute  the  best  estimate  of  the  TEC  along  the  line  of  sight 
to  the  satellite.  For  example. 


TEC(  vertical  )  =  W* 


+  (]-W)* 
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r(90°-EL)l 
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where  W  is  a  linear  weighting  factor,  and  i  =  1,5. 


(7) 


Kalman  filters  are  discussed  in  detail  by  Jazwinski  [11].  The  specific  filter  used  in  our  GPS 
ion  model  is  described  in  Appendix  B.  The  coefficients  o( ,  i  =  0,...,5  (defined  as  the  filter  state)  are 
estimated  in  the  Kalman  filtering  process.  The  observations  that  map  into  the  estimation  of  the  filter 
state  are  the  individual  TEC  measurements  and  their  corresponding  azimuths  and  elevations.  One 
set  of  measurements  exist  for  each  satellite  in  view  at  any  given  time.  For  example,  if  four  satellites 
were  viewed  at  time  fQ,  then  there  will  be  four  observations  of  the  TEC,  plus  their  associated 
azimuths  and  elevations,  that  will  be  used  to  compute  that  particular  filter  state  along  with  all  past 
measurements  of  TEC  that  have  been  appropriately  age-weighted.  The  TEC  measurements  used 
are  the  combined  pseudorange  and  integrated  phase  measurements. 


Figures  4(a),  (b),  and  (c)  show  a  comparison  of  the  vertical  TEC  values  above  Millstone  as  a 
function  of  the  time  of  day.  Figure  4<^)  shows  the  vertical  TEC  predicted  by  our  real-time  GPS 
model.  Figure  4(b)  shows  the  TEC  values  predicted  by  the  Bent  ionospheric  model.  Figure  4(c) 
shows  the  TEC  values  predicted  by  the  Slab  model  of  the  ionosphere  [12],  combined  with  a  real¬ 
time  measurement  of  the  foF2.  The  Slab  model  of  the  ionosphere,  combined  with  the  real-time  foF2 
information,  is  the  “old”  ion  model  used  in  the  real-time  satellite-tracking  program  at  Millstone.  As 
can  be  seen  in  these  three  figures,  on  average,  the  zenith  values  agree  among  the  different  models. 
The  diurnal  trend  of  the  TEC  is  clearly  evident,  although  the  GPS  model  shows  a  slightly  different 
behavior  than  the  Bent  and  Slab  models  predict. 
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5.  LAGEOS  TRACKING 


The  Millstone  radar  tracks  the  Lageos  satellite  on  a  regular  basis,  nearly  once  a  day.  Lageos 
is  a  satellite  with  a  perigee  of  approximately  5833  km  and  an  eccentricity  of  0.0047.  Lageos  is 
equipped  with  cube  comer  reflectors  and  thus  can  be  tracked  by  NASA’s  laser  ranging  system. 
These  lasers  have  a  ranging  accuracy  of  10  cm  or  better. 

The  Lageos  laser  ranging  data  are  used  to  monitor  the  calibration  of  the  Millstone  radar.  Using 
an  in-house  orbit  determination  program  called  Dynamo  [  1 3]  and  the  precision  laser  ranging  data, 
the  position  of  Lageos  can  be  predicted  to  within  20  cm.  By  comparing  our  measured  observations 
of  Lageos  with  the  predictions  of  Lageos,  an  estimate  of  the  system  bias  in  our  radar  can  be 
determined.  This  procedure  is  followed  at  Millstone  on  a  weekly  basis,  as  radar  calibration  is  an 
ongoing  process. 

When  the  real-time  GPS  ionospheric  monitoring  system  was  integrated  into  the  satellite¬ 
tracking  system  last  March,  we  frequently  switched  between  the  two  methods  of  estimating  the 
ionospheric  correction:  the  new  method  that  uses  the  GPS  Kalman  filter  and  the  old  method  that 
uses  real-time  foF2  information  combined  with  a  Slab  model  of  the  ionosphere  [  1 2].  We  then 
identified  which  Lageos  passes  used  the  new  GPS  ionospheric  model  and  which  used  the  previous 
ionospheric  model  to  determine  the  ionospheric  delay  term.  Figures  5(a)  and  (b)  show  the  difference 
in  the  estimates  of  the  average  Lageos  range  bias  for  the  two  ionospheric  models  and  for  several  of 
these  passes. 

Notice  the  difference  in  the  size  of  the  standard  deviation  between  the  data  using  the  GPS  ion 
model  (5.94  TEC  units)  versus  the  data  using  the  Slab  model  (22.02  TEC  units).  Clearly,  the  GPS 
ion  model  represents  a  significant  improvement  in  our  ionospheric  correction  estimate.  Interest¬ 
ingly,  the  residuals  that  use  the  old  ion  model  in  Figure  5(b)  seem  to  improve  during  the  latter  half 
of  the  year.  This  improvement  occurs  partially  because  fewer  passes  were  taken  using  the  old  ion 
model  and  partially  because  several  of  the  passes  contain  at  least  some  GPS  data.  To  qualify  as  being 
a  GPS  pass,  more  than  50%  of  the  data  in  a  pass  had  to  have  been  taken  using  the  GPS  data  to  correct 
for  the  ionospheric  term. 

It  is  also  worth  observing  that  several  of  the  Lageos  passes  were  taken  during  geomagnetically 
disturbed  time  periods  (defined  here  as  having  a  daily  Ap  greater  than  20).  In  fact,  half  of  the  Lageos 
passes  taken  using  the  GPS  ion  model  qualify  as  “disturbed”  time  periods.  It  is  extremely 
noteworthy  that  there  is  no  evidence  in  the  residuals  of  which  days  were  geomagnetically  disturbed. 

One  last  figure  is  presented  that  shows  the  dramatic  difference  in  the  prediction  of  the 
ionospheric  range  correction  between  the  GPS  model  and  the  old  ion  model  used  at  Millstone.  The 
data  were  taken  on  20  July  1 99 1 ,  a  day  when  the  daily  Ap  was  32  and,  in  fact,  the  Ap  had  been  elevated 
for  the  entire  prior  week.  In  Figure  6,  the  ionospheric  range  correction  estimated  by  both  the  GPS 
model  and  the  Slab  ion  model  is  shown  over  an  entire  pass  of  Lageos.  The  GPS  model  predicts  a 
range  correction  that  is  noticeably  less  than  the  old  ion  model,  particularly  at  the  beginning  of  the 
pass  that  corresponds  to  low  elevation  data.  The  Lageos  range  residual  that  was  calculated  using 
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the  GPS  estimate  of  the  ionospheric  correction,  is  the  bottom  line  on  the  axis.  The  Lageos  range 
residual  remains  close  to  zero  during  the  entire  pass.  This  residual  clearly  indicates  that  the  GPS 
estimate  of  the  range  correction  is  correct  and  is  a  significantly  better  estimate  of  the  true  ionospheric 
correction  than  that  of  the  previous  model. 
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FigureS.  (a)  Lageos  residual  using  GPS  ion  model.  Standard  deviation  =  5.94  TEC  units,  (h)  Lageos  residual  using  old 
ion  model.  Standard  deviation  =  22.02  TEC  units. 
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Figure  6.  Lageos  range  bias  versus  refraction  corrections:  GPS  and  old  ion  model  estimates,  20  July  1991 . 
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6.  CONCLUSION 


A  real-time  ionospheric  monitoring  system  has  been  developed  and  integrated  into  the 
Millstone  Hill  satellite-tracking  operation.  Initial  results  clearly  indicate  significant  improvement 
in  the  estimation  of  the  ionospheric  delay  term. 
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APPENDIX  A 

IONOSPHERIC  MAPPING  FUNCTION 


The  mapping  function  is  used  to  transform  the  zenith  or  vertical  path  delay  or  TEC  to  an 
arbitrary  zenith  distance.  The  model  assumes  a  spherically  stratified  ionosphere  with  a  simple 
height  distribution  of  electrons.  The  model  is  an  extension  of  the  Slab  model,  which  is  often  used 
to  interpret  and  apply  ionosonde  foF2  measurements  to  TEC  estimates.  The  distribution  of  electrons 
is 
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^max *s  l^e  •ower  boundary  of  the  slab, 

S  is  the  thickness  of  the  slab, 

^upper 's  toP  ionosphere, 

A bo„oni 's  bottom  of  the  ionosphere,  and 

A  and  B  describe  the  decay  of  the  electron  density  from  the  center  of  the  slab  to  the  “bottom”  and 
the  “top”  of  the  ionosphere,  respectively.  Nominal  values  can  be  Amax  =  300  km,  S  =  200  km,  /tupper 
=  1000  km,  ^bottom  =  100  km,  A  =  30  km,  and  B=  100  km.  The  model  is  normalized  to  a  zenith  TEC 
of  1 ,  and  the  mapping  function  can  be  applied  directly  to  calculate  the  TEC  along  a  path  at  a  given 
zenith  distance  and  range.  It  can  therefore  be  used  to  calculate  the  TEC  or  path  delay  to  a  point  in 
the  ionosphere. 

The  mapping  function  is  formulated  in  terms  of  the  zenith  distance  z,  the  range  to  the  satellite 
from  the  observing  station  r,  and  the  geocentric  height  of  the  station  r  .  The  height  of  the  target  h 
is  calculated  with 


h  =  Jr2  +  r2  +  2r  rcos(z)  -  r 

Vo  o  '  '  o 

Now  the  contribution  from  the  slab  is 


(A-2) 
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If  the  target  height,  /jmax  <  h  <  hmax  +  S,  then  this  contribution  is  evaluated  with  S  =  h  -  h  If  h 
<  hmax ,  there  is  no  contribution  from  the  slab. 

The  decay  of  electron  charge  from  the  slab  depends  on  the  parameters  A  and  B.  The  integrals 
can  be  done  analytically  and  are  taken  from  Gradshteyn  &  Ryzhik  [14],  Equations  2.275.4,  and 
2.266.  They  are 
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The  contribution  from  the  lower  ionosphere  is 
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The  contribution  from  the  upper  ionosphere  is 
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The  limits  of  integration  are  adjusted  in  the  obvious  way  for  an  object  within  the  ionosphere. 
Finally,  the  mapp.ng  function  is  normalized  to  the  zenith  using 
Z..(z,r)  +  Z  (z,r)  +  Z.  (r,r) 

slab'  ’  upper'  '  lower'  ’ 
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APPENDIX  B 

SEQUENTIAL  ESTIMATOR  OF  THE  CIRCUS  TENT  MODEL 


The  Circus  Tent  model  of  the  ionosphere  was  designed  specificall''  as  a  local  model  and 
follows  traditional  modeling  by  employing  a  mapping  function.  The  parameterization  also  provides 
a  comparison  with  other  data  (e.g.,  zenith  path  delay  from  an  ionosonde).  Six  disposable  parameters 
were  chosen  to  correspond  to  the  data  available  (i.e.,  a  four-channel  TI4100  GPS  receiver). 


The  Circus  Tent  model  represents  the  vertical  total  electron  content  (TEC)  as  a  function  of 
zenith  distance  (r  =  l/27t-ele)  and  azimuth  (Az).  The  half  space  centered  at  the  station  is  divided 
into  five  equal  pie-shaped  sectors  (quintets).  The  azimuth  of  the  five  sides  are  346°,  58°,  1 30°,  202°, 
and  274°.  These  sectors  orient  the  Circus  Tent  model  along  the  natural  magnetic  coordinates  of 
Millstone  Hill.  The  TEC  along  a  sector  edge  at  Azt  is  modeled  as 


~ 

+  a, 

- 

U/2-1 

(B  - 1) 


The  TEC  is  never  allowed  to  be  negative.  The  TEC  at  an  arbitrary  azimuth  is  obtained  by  linear 
interpolation  as 


TEC(r,  A:)  =  TEcJz,  /4ry)  +  TEC(z,  A:.)- TEc|r,  Az.} 


Az  -  Az . 
_ J 
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where  Azj  <  Az  <  Azr  Note  that  aQ  corresponds  to  the  zenith  TEC.  The  quadratic  dependence  of 
zenith  distance  was  found  to  be  a  better  model  than  a  linear  dependence.  This  is  probably  because 
the  mapping  function  is  very  nearly  linear  and  already  models  the  ionosphere  quite  well. 


Three  separate  steps  are  involved  in  the  sequential  estimate  of  the  six  model  parameters  that 
define  the  Kalman  filter  state.  The  first  step  is  correcting  possible  cycle  slips  in  the  carrier  phase 
data,  which  is  in  addition  to  the  cycle  slip  corrections  determined  in  the  preprocessing,  as  described 
in  Section  3.  The  second  step  is  using  the  carrier  phase  to  smooth  the  pseudorange.  Finally,  the  third 
step  is  the  prediction  of  the  state  at  the  observation  time  and  the  correction  of  the  state  using  the 
observation  vector. 


CYCLE  SLIP  CORRECTION 


The  cycle  slip  correction  models  the  carrier  phase  <pL]  and  <pL 2  as  rational  polynomials 


<P(X) 
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+  a^x  +  a^x 
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(B  -  3) 


The  existence  of  a  cycle  slip  in  L 1  and  L2  is  tested  as  follows.  First,  the  polynomial  coefficients  a0, 
a,,  a2,  and  o3  are  computed  by  least  squares  from  the  previous  20  measurements  of  carrier  phase. 
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The  first  20  observations  in  a  sequence  are  assumed  to  be  error  free.  The  independent  variable  of 
the  fit  is  centered  and  scaled  to  the  data  interval  with 


x.  = 
/ 


i  ^mid 


max  min 


(B  -  4) 


where  tj  is  the  observation  time,  tmid  is  the  average  of  all  the  times  in  the  interval,  rmin  is  the  first  time 
in  the  interval,  and  /max  is  the  last  time  in  the  interval.  The  LI  and  L2  phases  are  then  predicted  for 
the  next  observation.  The  difference  in  predicted  and  observed  phase  for  LI  and  L2  are  converted 
to  cycles  using  0. 1 9029  m  and  0.2442 1  m  as  the  wavelengths.  Then  the  difference  in  predicted  phase 
is  compared  with  the  difference  in  observed  phase.  If  the  difference  in  phase  difference  exceeds 
4.5  cm  and  integral  cycle  slips  are  predicted,  a  cycle  slip  is  declared.  The  integral  number  of  cycles, 
expressed  in  meters,  is  then  added  to  all  subsequent  observations  of  LI  and  L2  as  appropriate.  The 
cycle  slip  correction  is  cumulative.  However,  if  three  consecutive  observations  have  a  cycle  slip  or 
the  number  of  cycles  slipped  exceeds  100,  the  process  aborts.  The  phase  accumulation  and 
smoothing  of  pseudorange  data  for  the  particular  satellite  are  then  restarted. 


PHASE  SMOOTHED  PSEUDORANGES 


The  difference  in  the  observed  pseudorange  A r,  corrected  for  satellite  and  receiver  bias  (de¬ 
fined  in  Section  2),  is  a  direct  measurement  of  the  TEC  along  the  LOS  to  the  satellite.  The  differences 
in  the  observed,  cycle  slip  corrected,  carrier  phase,  and  A0  (in  meters)  are  different  from  the  LOS 
TEC  by  a  constant  (i.e.,  the  LOS  TEC  at  the  beginning  of  the  set  of  data).  Therefore,  the  smoothed 
data  A<f>  for  the  nth  data  point  is  taken  as 
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Some  data  screening  is  possible.  First,  a  pseudorange  difference  or  carrier  phase  difference  that  is 
negative  or  more  than  50  m  is  declared  an  error,  and  the  observation  is  skipped.  Five  consecutive 
errors  cause  the  phase  accumulation  and  smoothing  of  pseudorange  data  for  a  particular  satellite  to 
be  restarted.  In  addition,  for  use  in  the  Kalman  sequential  estimator,  each  observation  has  an 
associated  uncertainty  based  on  the  standard  error  of  fit  and  an  observation  noise  of  0.20  m. 
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SEQUENTIAL  ESTIMATOR 

The  sequential  estimator  developed  here  follows  Jazwinski’s  theoretical  framework  [1  1).  The 
sequential  estimator  processes  the  smoothed  pseudorange  data,  corrected  for  satellite  bias,  as 
measurements  of  LOS  TEC.  The  measurement  y  is  a  vector  of  up  to  four  contemporaneous  Ap's 
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from  satellites  1 , 2,  3,  and  4.  The  data  are  used  to  estimate  the  ionospheric  state  x  made  up  of  the 
six  Circus  Tent  parameters 
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where  aQ  is  th  z^"1  \  TEC,  and  ar  i  =  1,  2,.. .,5  are  the  slopes  along  each  of  the  quintet  sector 
boundaries.  2  ,ie  system  model  is 


J>  x  .  +  T  .w  .  . 
7  7  77+1 


(B  -  9) 


The  sta»'  transition  matrix  <I> .  =  <t>  (/  tj)  that  predicts  the  slate  at  r+|,  given  the  state  at  /,  is 


'10  0  0  0  0' 

0  /  0  0  0  0 

0  0  /  0  0  0 

0  0  0  /  0  0 

0  0  0  0  /0 

0  0  0  0  0  / 
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x  .  .  =  4>  .x  . 
7+1  7  7 


7+1  7 


(B  -  1 1) 


Therefore,  the  zenith  TEC  is  constant,  and  /  =  e  ^  ’  >s  chosen  so  that  the  TEC  slope  becomes 

zero  after  some  time.  These  are  the  best  estimates,  in  the  absence  of  data,  to  correct  the  state  estimate. 
x  is  chosen  to  be  180  min. 
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The  process  noise  vt^  represents  ionospheric  model  uncertainties,  and  dynamic  and  stochastic 
ionospheric  effects  are  excluded  from  the  model.  These  effects  are  assumed  to  be  a  Gauss-Markov 
process  with  zero  mean  and  a  variance  qr  Now,  each  element  of  the  state  x-  is  assumed  to  have  an 
uncorrelated  and  independent  Gauss-Markov  process.  Therefore,  T  is  the  identity  matrix,  I  (T  =  I). 

The  predicted  state  error  covariance  P  depends  on  0,  T,  and  Q  where  Q  is  the  diagonal  matrix 


0  0  0  0  0 

0  <7,  0  0  0  0 

0  0  q2  0  0  0 

0  0  0  q  0  0 

0  0  0  0  q4  0 

0  0  0  0  0  q5 


(B  -  12) 


The  ith  element  qi  is  the  variance  of  the  Gauss-Markov  process  of  the  /th  parameter  over  the  yth  interval. 
q,  is  nominally  set  to  10"6.  So  the  predicted  error  covariance  is 


p  =<p  P0T  +  r q  trT  . 

7+1  777  77+1  7 

The  measurement  update  uses  the  measurement  matrix  of  partial  derivatives: 
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where  = 
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the  measurement  error  covariance  R  is  taken  from  the  computed,  smoothed  pseudorange  error 
estimates  cx  derived  in  Equation  (B-6)  as 


CTj  0  0  0 

0  ct2  0  0 

0  0  cx3  0 

0  0  0  a . 

4 


(B-  15) 
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The  measurement  update  uses  a  fading  memory  to  allow  real  changes  in  the  ionosphere  to  update 
the  model  and  prevent  the  consequent  filter  divergence.  The  fading  memory  is  implemented  by 
modifying  the  error  covariance  by  the  factor /,  as  seen  in  Equation  (B- 1 0).  The  fading  memory  time 
constant  r  is  assumed  to  be  the  same  as  ionosphere  decay  time  (180  min),  although  there  is  no 
physical  basis  for  this  assignment.  With  this  formalism,  the  Kalman  gain  is 


K 

j 
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Then  the  corrected  state  estimate  is 


=  x  ,K 

7+1  7 


y  -  M  .x  .  . 
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and  its  associated  covariance  is 


(B  -  17) 
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This  sequential  estimator  has  been  developed  to  be  data  driven  and  to  make  few  physical 
assumptions  about  the  ionosphere.  It  works  well  under  a  variety  of  different  conditions.  For 
example,  this  estimator  works  well  with  one  or  more  GPS  satellites  in  view,  and  it  continues  to  work 
when  the  constellation  of  GPS  satellites  in  view  changes.  This  estimator  also  seems  fairly  robust 
during  disturbed  ionospheric  conditions.  Finally,  valid  recovery  of  ionosphere  parameters  was  seen 
even  in  the  presence  of  selective  availability  (S  A),  although  there  are  clearly  carrier  dither  sequences 
that  can  defeat  the  cycle  slip  detection  and  correction.  This  method  of  estimation  can  easily  be 
expanded  to  make  use  of  the  following:  more  parameters,  more  observations  to  update  the  state,  the 
inclusion  of  a  physics-based  model  for  the  state  transition  matrix,  and  better  models  for  the 
measurement  error  covariances  and  the  process  noise. 
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